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Superconvergence and recovery type a posteriori 
error estimation for hybrid stress finite element 

method* * 

Yanhong Baij Yongke Wuf Xiaoping Xie§ 


Abstract 

Superconvergence and a posteriori error estimators of recovery type 
are analyzed for the 4-node hybrid stress quadrilateral finite element 
method proposed by Plan and Sumihara (Int. J. Numer. Meth. En- 
grg., 1984, 20: 1685-1695) for linear elasticity problems. Uniform 
superconvergence of order with respect to the Lame 

constant A is established for both the recovered gradients of the dis¬ 
placement vector and the stress tensor under a mesh assumption, where 
a > 0 is a parameter characterizing the distortion of meshes from par¬ 
allelograms to quadrilaterals. A posteriori error estimators based on 
the recovered quantities are shown to be asymptotically exact. Nu¬ 
merical experiments confirm the theoretical results. 

Keywords: linear elasticity, hybrid stress finite element, superconver- 
gence, recovery, a posteriori error estimator 


1 Introduction 

Assumed stress hybrid finite element method (also called hybrid stress 
method) pioneered by Pian [30] is known to be an efficient approach in 
the analysis of elasticity problems (cf. [HU [321 IMl IMl E51 ESI SSI ESI ES] ) • 
One main advantage of the hybrid method lies in that, the method allows 
for piecewise-independent approximation to the stress solution and, through 
local elimination of the stress unknowns, finally leads to a symmetric and 
positive definite discrete system of unknowns of displacements. In [33| Pian 
and Sumihara derived a robust 4-node hybrid stress quadrilateral element 
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(abbr. PS) through a rational choice of stress terms, where the continuous 
piecewise isoparametric bilinear interpolation is used for the displacement 
approximation. We refer to [48] for the analysis of uniform convergence and 
a posteriori error estimation for the hybrid stress quadrilateral elements 
proposed in [33l US] . 

As an active research topic, super convergence of finite element solutions 
to partial differential equations has been studied intensively for conforming, 
nonconforming and mixed finite element methods (see, e.g., books mm 

EH EH SH Sni EZ] and papers [H El EH EH ESI ESI ESI E3 ESI SSI SSI 

SH ESI El ES]). Based on theory of superconvergence, a posteriori error 
estimation of recovery type has attracted more and more research interests 
in recent two decades. The most representative recovery type error estimator 
is the Zienkiewicz-Zhu (ZZ) estimator based on gradient patch recovery by 
local discrete least-squares fitting [58l [59] . The method is widely used in 
engineering practice for its robustness. Superconvergence properties of the 
ZZ patch recovery were shown in [5HE9| for rectangular and strongly regular 
triangular meshes, respectively. The work of mm introduced a recovery type 
error estimator based on global L^-projection with smoothing iteration of 
the multigrid method, and established asymptotic exactness in the H^-novm 
for linear element under shape regular triangulation. By using the result in 
[6], a new theoretical justification was given in [l6| for the ZZ estimator. 
A polynomial preserving gradient recovery (PPR) method was proposed in 
[531154] which is different from the ZZ gradient patch recovery method [58] . 
In [IH some patch recovery methods were proposed and analyzed for finite 
element approximation of elasticity problems using quadrilateral meshes. 

So far, to the authors’ knowledge, there is no super convergence analy¬ 
sis for the hybrid stress finite element method for the elasticity problems. 
This paper is to establish superconvergence for the Pian and Sumihara’s 
hybrid stress quadrilateral element [33] • We shall derive the uniform su¬ 
perconvergence with respect to the Lame constant A for both the recovered 
displacement gradients and the recovered stress tensor, and show that the a 
posteriori error estimators based on the recovered quantities are asymptot¬ 
ically exact. 

The rest of the paper is organized as follows. Section 2 introduces the 
model problem and its weak form. Section 3 shows the hybrid stress finite 
element discretization and some preliminary results. Section 4 analyzes the 
superconvergence of the hybrid stress method. Section 5 is devoted to the 
recovery of the displacement gradients and the stress tensor, as well as the a 
posteriori estimation of recovered type. Finally, Section 6 provides numerical 
results. 
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2 Model problem 

Let 0 C be a bounded polygonal domain with boundary dQ. We con¬ 
sider the following linear elasticity problem with homogeneous displacement 
boundary condition: 


—diver 

a 

u 


f in 0, 

Ce(u) in 0, 

0 on r := dfl, 


( 2 . 1 ) 


where Q C is a bounded polygonal domain, a € denotes the sym¬ 

metric stress tensor field, u € the displacement hied, e(u) = ^ (Vu -|- (Vu)^) 
the strain tensor, f € the body loading density, and C the elasticity mod¬ 
ule tensor with 

Ce(u) = 2/re(u) -|- AdivuX. 

Here Z is the 2x2 identity tensor, tr((T) the trace of the stress tensor a, 
and /r, A the Lame parameters. 

We introduce some notations as follows. For an arbitrary open set T, 
we denote by H^{T) the usual Sobolev space consisting of functions defined 
on T with derivatives of order up to k being square-integrable, with norm 
II • Ilfc^'T and semi-norm | • \k,T- In particular, H^(T) = L?‘{T). When T = Q, 
we abbreviate || • ||a;,q and | • \k,n to || • ||fc and | • \k, respectively, and denote 
II ■ II II ■ llo- We use the same notations of norms and semi-norms as above 
for corresponding vector or tensor spaces. For any vector a = € M"", 


we denote 


and 


\a\ 


loo := max 

l<i<n 


ai 


,I«IIf := E 

Vj=i 

Throughout the paper, we use notation a < 6 (or a > 6) to represent that 
there exists a constant C, independent of mesh size h and the Lame constant 
A, such that a < Cb (or a > Cb), and use a ^ b to denote a < b < a. 

Dehne the spaces 


S := <^ r G 


2x2 

sym 


), [ tr(r) 
Jn 


= 0 


V := = {v € : v|r = 0}, 

where L^(H; denotes the space of square-integrable symmetric tensors, 

and tr(r) := th -I-T 22 the trace of tensor r. Then we have the following weak 
problem for the system m- Find (u, u) G S X V such that 

f a((T,r) -b 6(t,u) = 0 for all r G 51, , . 

\ b{a,v) = F{v) for all v G V, ^ ' 
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where 


6(t,v) = -/' T:e(v), F(v) = - / f • v. 

Jo, J n 

It is well-known that the weak problem (12.211 admits a unique solution. 


3 Hybrid stress finite element discretization 


3.1 Geometric properties of quadrilateral meshes 

Let {Th}h>o be a- partition of O by convex quadrilaterals with the mesh 

size h := max hx, where hx is the diameter of quadrilateral K G Th- 
KdTh 

Let Zi{xf ,yf) and Zi{^i,'qi) for 1 < z < 4 be the vertices of K and the 
reference element K = [—1,1]^ (cf. Figure [3d]), respectively. There exits a 
unique invertible bilinear mapping Fx : K ^ K that maps K onto K with 
Fx{Zi) = Zi. The mapping Fx is of the form 




an -h af ^ -h a^r] -h a^ 2 ^ri \ 

+ bfi + b^'n + bUri ) 


(3.3) 


where zy € [— 1 , 1 ] are the local coordinates and 


/ 


b^ 

\ 


/ 

1 

1 

1 



/ 

Xi 

Vi 

\ 


of 

6f 


1 


-1 

1 

1 

-1 



X2 

2 /f 



of 

b^ 


“ 4 


-1 

-1 

1 

1 




2 /f 


V 

of 

6f 

/ 


V 

1 

-1 

1 

-1 y 


V 

d/4 

vf 

/ 


In the following we may omit the superscript K of the above notations if 
there is no confusing. 

The Jacobi matrix and Jacobian of Fx are respectively given by 


DFx{i,y) 


dx 

dx 


dy 

dy 

dy 

di 

dy 


ai + ai2'n 

bi + bi2r] 62 bi2C 


(3.4) 


Jk{C, v) = det{DFx) = Jo + JiC + «l 2 z?, (3.5) 

where 


Jo — aib2 — a2bi, Ji — 01612 — 01261 , J 2 — 01262 — 02612 . 


It is easy to obtain the inverse of the Jacobi DFx with 


DF-^oFx{i,y) 


di di 
dx dy 
dy dy 
dx dy 


1 f b2 + 612? —O2 — O12C \ 

Jk V ~bi - buv Oi -k oi2zy J ' 

(3.6) 
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Figure 3.1: Bilinear transformation Fk maps reference element K (in the 
left) to element K (in the right). 


Throughout this paper we assume the partition Th is shape regular in 
the following sense [50] : There exist a constant g > 2, independent of h, 
such that for all iF € 7/i it holds 


hx < QPk- (3.7) 

Here px ■= min Pi, with pi being the diameter of the largest circle inscribed 
l<j<4 

in Tj, the sub-triangle of K with vertices Zj_i, Zi and Zj+i (the index on 
Zi is modulo 4) for i = 1, • • • ,4. 

We introduce several additional mesh conditions which will be used in 
the forthcoming analysis of super convergence (Section S]). 

• (MCI) Diagonal condition: There exists a constant a > 0 such 
that for any quadrilateral K G Th, the distance, dx {dx = IO1O2I = 

+ 6 ^ 2 )) between the midpoints of the diagonals of K (See Figure 
El]) satisfies 

dx = 0{h]+'^). (3.8) 

• (MC2) Neighboring condition: For any two quadrilaterals ii'i, K 2 G 
Th sharing a common edge, it holds, for j = 1 , 2 , 

= af^{l + 0{h%^ + +h%^)). (3.9) 

Remark 3.1. Diagonal condition (MCI) is also called (1 + a)-section 
condition (cf. Note that K is a parallelogram if and only if dx = 0, 

which means a = -|-oo. When a = 1, (MCI) is the Bi-Section Condition or 
condition B 

Remark 3.2. Th is said to satisfy Jamet condition m if there exists a 
constant r > 0 such that hx < rpx holds for any quadrilateral K £ Th, 


5 












where px is the diameter of the largest cirele inseribed in K. As shown in 
]23^ . if both Jamet condition and Diagonal condition (MCI) hold, then 
Th is shape regular for suffieiently smal h. 

In view of the shape regularity condition (j3.7|l . it is easy to obtain the 
following estimates for the Jacobian Jx given in (13.51) . 

Lemma 3.3. For any K ^Th it holds 

Jft- PS Jo ~ hj^. (3.10) 

Further more, if Diagonal condition (MCI) is satisfied, then it holds 

max{|Ji|, IJ 2 I} pa/i^". (3.11) 

3.2 Pian-Sumihara’s hybrid stress finite element method 

In view of the mapping Fx, for any function w{^,rj) on K we define 
function w{x,y) on K gTh with 

w{x, y) := w(^, rj) or equivalently w := w o Ff)^. 

In Pian-Sumihara’s hybrid stress finite element (abbr. PS element) 
method [33] for the problem (12.2p . continuous piecewise isoparametric bi¬ 
linear interpolation is used for the approximation of displacement, namely 
the displacement approximation space V/j C V is taken as 

V;, := Sh X Sh 


with 


Sh = {v e : V = v\x o Fx e span{l,r/,^r/}, for all K G %}■ 

To describe the stress approximation of PS element, we abbreviate the 


symmetric tensor r = 


Til Ti2 


T12 T 22 

of PS element is of the following form on K: 


to T = {tii,T 22 ,ti 2 )'^■ The stress mode 


T = 


Jll \ 

f 1 

0 

0 

7 ] 

\ 


T22 = 

0 

1 

0 

bl 


=; d/3^, /3^ G (3.12) 

ri2 / 

VO 

0 

1 

ai 1 

62^ / 



Then the corresponding stress approximation space, C S, for PS element 
is given by 


S/j := {r G S : f = t\x ° Fx is of the form (I3.12p for all K G %,}■ (3.13) 
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As a result, the PS element method for the problem (j2.2ll is given as follows. 
Find {ah, Uh) G S/j x V/j such that 

f a{ah,T) + b{T,Uh) = 0 for all r G Sh, 

\ b{ah,v) = F(v) for all V G V/j. ^ ‘ 

Let (<T, u) G P| S) x (V be the solution of the 

problem (I2.2I1 . It has been shown in [38] that the following uniform error 
estimate holds for the PS element method: 

||cT - + |u - Uh\i < h(||u ||2 + Iklli). (3.15) 


4 Super convergence analysis 


4.1 Preliminary results 


We recall v{x, y) := v{^, r/) = v o Fk {x,y). Some calculations show 


d^v 

d'^v 

d^dy 


( d ^ \ F 

{ai + ai 2 r])-^ +ibi+bi 2 r])-^] v, r = 1,2, (4.16) 

dv dvf d , \ 9 \ 

““s + + (<“■ + 

X ((02+ ai2^}^ + (f.2 + 6 i 2^)^) If. (4-17) 


In light of these two relations and Lemma f3.3[ we easily derive the following 
lemma. 


Lemma 4.1. For all K ^Fh and v G H^{K), it holds 

dv 


dv 

+ 

d»4- 

o,k 


d'^v 


5^2 

o 

+ 


dr] 

dH 


^ \v\l,K, 


drf 


0,K 


0,K 


< 


hK\v\ 


2,K- 


(4.18) 

(4.19) 


In particular, if Th satisfies Diagonal condition (MCI), then it holds 
di^v 


d^drj 


^ hK\''^\l,K + hK\v\2,K- 


(4.20) 


0,K 


Let G V/i be the piecewise isoparametric bilinear interpolation of 
u G V P|(i4^(n))^, then it holds the following estimate: 

||u - u^llo.E-+ hii-|u-< /i^||u||2,2;v, for all AT G 7^. (4.21) 

Let a^ G T,h be the projection of a G S in the a{-, •)-inner product, 
namely a^ satisfies 


a{a^, t) = a{a,T) for all r G S/j. (4.22) 
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Thanks to ()3.12p and ()4.22p . we obtain, for all K € Th, 


cjV = 


^ I C“V with Hk ■■= / ^ C"M, 

Jk Jk 


(4.23) 

(4.24) 


[ = 0 - 

Jk 

In addition, we have the following lemma. 

Lemma 4.2. Under Diagonal condition (MCI), for all K ^ Th it holds 

\W - o-^\\o,K < hK\W\\i,K, (4.25) 


(d - M) 


IK 


<h 


kWTTk- 


(4.26) 


Proof. Let € T^h be the projection of a with 


d'-T=a-T, for all T G S/j. 

In Jn 


Then we have 


cr — cr 


ill < 




and 


a\K = / Al a with Hk '■= / Al A 


JK Jk 

for sl\ K ^ Th- By triangle inequality, it holds 

||<i — (T^W < ||<T — d’^ll + lld-^ — a^W < h||fT||i + — cT^I 


(4.27) 


(4.28) 


(4.29) 


We turn to estimate \\d^ “ = {^k&Th 

(I4.22P and (|4.28l) . some calculations yield 


\ 1/2 

|d^ —. In view of 


H-^A^C-^ = — 
^ 4Jo 


/I 0 0 \ 

0 10 
0 0 1 

dAiV dA2'n disg 

V d^2f. d^sf, ) 


+ h.o.t., 


(4.30) 














- 4Jo 


1 

0 

0 


3v 


4^ 

“2 I “2 


0 

1 

0 


4^ 

°1 

~ 6 ? bF 


3S 


0 

0 

1 


3-^j? 

a,^ a^ 




2 

^—:t 


+ h.o.t., (4.31) 


where 
dii = 

f^5i = 


1 + 1 

3 ( £2 


f^42 — 


2/i+A 

1 + ^ 


2 ’ 


1 + ^ 


2\2 ’ 


d^2 — 


3 1- 


“i A 

2/i+A 


1 + ^ 


2 \ 2 


IT 


f^43 — 


1^53 — 


12 (^ + A) 


ai 


2 , + A 

12(m + A) e 


2 ’ 


2/x + A 


1 + ^ 


2 \ 2 ’ 


and in each of the above two relations h.o.t denotes a different higher-order- 
term matrix of the form 

h.o.t = ^(hjj(C,r?))5x3 with max max \hij\ < h%. (4.32) 

Jo 


Obviously, it holds 

max \dij \ « 1. (4.33) 

4<i<5, l<i<3 

Denote QkO' := then a combination of (13.101) . (13.111) . (14.221) 

and ()4.28p - ()4.32l) leads to 


\a^ — d^l 


,K = \\A -H],^A^')a\\o,K 


hK\ 




K 


hK\ 

Ik 

^ dh\\i,K, 


H-^A^C-^ - H-^aA {a - Qko)\ 


which, together with (|4.29l) . indicates the desired result (|4.25|) . 
The thing left is to prove (I4.26p . From (|4.24l) it follows 


0 = 


IK 


[a -a^) = [ Jk{^ - o-^) 

Jk 


= Jo ia-a^) + .Ji C{d-a^) + J2 r]{a-a^), 

Jk Jk Jk 
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which, together with (|3.ini) - ()3.1ip and ()4.25D . implies 



cr^ 


< 


< 

r\j 

< 

< 

rsj 


h%\\a — 


+ 


y / ??(d - 
Jo Jk 


^0,K 


hK\\Al,K- 


hx^xWcr - fT* 


lo,;^ 


□ 


For any K & Th, follow m to define the modified partial derivatives 
1 ^, and the modified strain tensor e(v) as 


{Jk-^\k o FK){i,r]) 
{JK^WoFKmv) 


(9y(0,0) dv 9y(0,0) dv dv dv 

drj dr] ^9^ ^ dr]' 

dx{0,0) dv dx{0,0) dv _ dv dv 

dr] d^ d^ dr] d^ ^ dr]' 


e(v) 


dvi 

dx 


Ifdvi I 8v 2\ 
\ 2^ dy ' dx r 


1 f dvi I 
2\ dy dx ) 


dv2 

dy 


(4.34) 


respectively. By the definition of e(v) it is easy to derive the following result. 


Lemma 4.3. Under Diagonal condition (MCI), for all v G V/i and 
K £Th it holds 

||e(v) -e(v)||o,ii' < 

Define the bubble function space \\ as 

\\ := jv'’ G ■■ v^(?, g) = v^Ie- oFk ^ span{^^ - 1, - 1}^, for all iF G y| . 

Then it is easy to verify that the PS stress mode (13.1211 satisfies the relation 
(see [37]) 

f e(v^) • r = 0, for all G V^, r G K eTh- (4.35) 

Jk 


4.2 Superconvergence analysis 

Define two functions 

E{Q ■= (ie - 1), F{r,) := - 1). 

Obviously it holds 

E'iO=^, E''{C) = 1, E'{r])=g, F" (g) = 1. (4.36) 
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Lemma 4.4. Under Diagonal condition (MCI) and Neighboring con¬ 
dition (MC2), for any g G and v G Sh it hold 

< h{h‘^\g\i + h\g\ 2 )\v\i, (4.37) 

Z] ~ h{h'"\g\i + h\g\2)\v\i. (4.38) 

Proof. We only give the proof of the first inequality, since the proof of the 
second one is similar. For any K G Th, g G and v G Sh, by (|4.36H . 

integration by parts, Cauchy-Schwardz inequality and Lemma l4. 11 we have 



+ {0{h]^°)\g\i^K + 0{h'j^)\g\2,K) \v\i,K, 


where lu and li are the upper and lower edges of K (see Figure IHTTI) . If the 
edge lu C (90, then the second term of the last equality in (I4.39p vanishes 
due to the homogeneous Dirichlet boundary condition, i.e. v\dn = 0. If lu is 
an interior edge of the partition 7h, we assume lu is shared by two elements, 
K and K^, of Th- By Neighboring condition (MC2) we have 

-/ixj = 0(/i'+“), 

then, from trace inequality and inverse inequality, it follows 

\hK-hKj\lu\f E{({s))^^ds 

J lu 

< h^~^°'{\g\i^K + hK\g\2,K)\v\i,K- (4.40) 

The above arguments also apply to the edge k. As a result, a combination 
of ()4.39p - ()4.40p yields the desired estimate ()4.37p . 

Similarly we can obtain (j4.38p . □ 

Lemma 4.5. Under Diagonal condition (MCI) and Neighboring con¬ 
dition (MC2), for a G 77^(0, M^^^) D S and u G (H^(Q))^ ClVit holds 

a{a^ — (j,t) = 0, for all r G T,h, (4-41) 

b{a-a^,v) < (/i^+"||(t||i +/i^|cr| 2 ) |v|i, for allv gNh, (4.42) 

6 (r,u-u^) < /i^+"||u|| 3 ||r||, forallTGT.h. (4.43) 


(e,i)de 


(4.39) 
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Proof. The relation (|4.4ip follows from (|4.22p . i.e. the definition of . 

Now we prove the estimate (|4.42p . For any v € V/j, we decompose it 
as V = vi + V2 with vi = vi|i^ o G span{l, r/}^, V2 = V2|ii- o Fk ^ 
span{,f7/}^, then it holds 



=: Ii + /2 + Is. (4.44) 


We note that Jft'e(vi) is a constant vector on K by the definition (I4.34p . 


'i'hns. in view of Lemmas 14.1114.21 we have 


\h\ = 

KGTh 

1 Y ■ [{a-a^)\ 

Ken 

< 

hKh\\l,K\\JKK^l)\\Q^K 



K&Th 


< 




K€Th 


< 


(4.45) 


For the term I 2 , from Lemmas 14.2114.31 it follows 

KeTh 


We turn to estimate I 3 . Denote V 2 = '^ 2 \k ° Fr ='■ and 

cr^ =: AI3^. Then, by (I4.23P and (|4.3Up . we have 


/ / 1 
0 


= Hr^[ 
Jk 


A^C-^a = 


IK 


1 

4^ 


V 


0 


0 

1 

0 


0 \ 
0 
1 


d4l'n d42'n d^sT] 

V d^if d^2f d^zi J 


\ 


+ h.o.t. 


cr. 
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Neighboring condition 

(MC2), Lemma [O] and (I4.32I) - (I4.33I) . yields 


which, together with uq = ^ 7 ^, 

j 

cr^ •e(v2)| = I (/3^)^ / ^(Jxe(v2)) 


E 

KeTh' 


K 


K&Th 


( 


E w ('»')' 




0 
0 
0 

, i>i 

ai 

h 


\ 


< + /i^|c 7|2) |v|i. 


(4.47) 


Similarly, since 

(T • e(v 2 ) = a- ■ (Jii-e(v 2 )) 

Jk 


JK 
it follows 


= 

Jk 


uohv - uoh^ 

-voa2r] + voai^ 

{vob2 - uoa2)r] + (liooi - vah)^ 


Y [ 0 -• e(v 2 )| < +/i^kb) |v|i, 

K^Th 


which, together with ()4.47D . yields 

\h\ < (/i^’^"|cr|i + /i^kb) |v|i. 


(4.48) 


As a result, the inequality (14.4211 follows from (I4.44I1 - (I4.46I1 and (I4.48p . 
The thing left is to prove the estimate (|4.43l) . Denote 


Xn : = 


b2 

-bi 

0 

0 ^ 

i / 

( bi2^ 

-buV 

0 

0 

0 

0 

-02 

Ol 

, Xi : = 

0 

0 

-012? 

012?/ 

-02 

Ol 

b2 

bi J 

1 \ 

V -012? 

012?/ 

612? 

bi2ri 


Vu : 

/ du 

du 
’ drj ’ 

dv dvY 
di' dr]) 

for u = 

{u,v), 




and let G be such that u-^ + u^ is the piecewise quadratic interpolation 
of u in the local coordinates rj. Thanks to the relation (14.3511 and the 
interpolation theory by [ 2 ] , for r G S/j it holds 


L 


e(u — u^) • T = 


IK 


< 


XoV(u — u'^ — u^) + AiV(u — u^)j ■ fd^dr] 

h - + h]+^\{i - Iloilo, J 


KK 


Then the desired inequality (I4.43P follows. 


□ 
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We are now in a position to state the following super convergence results 
for the hybrid stress method (j3.14l) . 

Theorem 4.6. Let (u, u) € f] S x V and {ah, € 

X V/j be the solutions of the problems (12.21) and (|3.14|) . respectively, and 
let € Yh be the isoparametric bilinear interpolation of u and a^ G S/j 
be the projection of a defined in (I4.22D . Then, under Diagonal condition 
(MCI) and Neighboring condition (MC2), it holds 

\Wh-cr^\\ < +/i^|(t|2 , (4.49) 

Iw-u'^li < /i^+"(||u||3 + ||cr||i) +/i^|(t| 2 . (4.50) 

Proof. From (j2.2l) and (j3.14p we easily obtain the error equations 

a(cr — cj/i, r) + 6(r, u — u/j) = 0 for all r G S/j, (4-51) 

b{a — ah,'v) = 0 for all v G V/i, (4-52) 


which, together with the discrete inf-sup condition for b{-,-) (cf. [H])! indi¬ 
cates 


Iw - u^li < sup 


b{T, VLh - u^) _ b{T, uh-u) + b{T, u - u^) 


= sup 


tg-eh Irll 

6 (r, u-u^) + a{a - ah, r) 

= sup -j—-j- 

tgEh Irll 

b{T, u-u^) + a{a - a^r) + a{a^ - ah, r) 

= sup -Tj—j-j- (4.53) 

TdEh Irll 


and 


Wh - < sup 


b{ah-a\v) b{a-a^,v) 


vGVh |v|i 


= sup 

vGVh |V|l 


(4.54) 


Then the desired estimates follows from the above two inequalities and 
Lemma 14.51 □ 


5 A posteriori error estimation of recovery type 


As shown in the estimate (j3.15D . the finite element solution {ah, u/^) of PS 
hybrid stress method (I3.14p is only of first order accuracy. We shall show in 
Subsections 15.1115.21 that, by using the recovery techniques of [2^ fiTl [45] . 
one can obtain recovered displacement gradients and stresses of improved 
accuracy, i.e. Then, in Subsection 15.31 we shall show the 

asymptotical exactness of the a posteriori error estimators based on the 
recovered quantities. 
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5.1 Gradient recovery by PPR 

We follow the polynomial preserving recovery method (PPR) proposed 
in [SHlEalll] to construct the recovered displacement gradients 

GhUh = {Ghui, Ghulf. (5.55) 

Here the gradient recovery operator Gh ■ ^ x Sh is dehned as follows 

[55] : Given function € S^, first define GhVh at all nodes (vertices) of the 
partition Th, and then obtain GhVh on the whole domain by interpolation 
using the original nodal shape functions of Sh¬ 
in PPR the values of GhVh at all vertices of Th are determined through 
the htting method. In fact, let Zi{xi,yi) be any interior vertex of Th, and 
let uji be a patch which consists of elements sharing the vertex Zi, i.e. 

LJi := ^ Th '■ Zi IS a vertex of K}. (5.56) 

For convenience all nodes on Qi (including Zi) are denoted by ^ij, j — 
1,2,-■■ ,n{n ^ 6). We use local coordinates (x,y) with Zi as the origin, 
i.e. (x,y) = where h := hi denotes the length of the longest 

element edge in the patch w,. The fitting polynomial is 

P2{x,y]Zi) = P^c (5.57) 


with 


P = {l,x,y,x^,xy,'f‘)'^, c = (ci,/iC 2 , hcs,/i^C 4 ,/i^cs,/i^ce)^. 


The coefficient vector c is determined by the linear system 

(5.58) 

fi \ 
yl 

yl / 

Finally, define 

GhVh{Zi) := Vp2(0, 0; Zi). (5.59) 

As shown in [Sam], under Diagonal condition (MCI) and Neigh¬ 
boring condition (MC2), the gradient recovery operator Gh is a bounded 
linear operator on the isoparametric bilinear displacement finite element 
space Vh = Sh X Sh T the followng sense: 


Q'^Qc = Q'^hh, 

where = {vh{Zii),VhiZi 2 ), • • • , Vh{Zin))'^ and 


/ 1 

Xl 

yi 

xl 

xiyi 

1 

X2 

m 

xj 

X2m 

V 1 

Xn 

ijn 

xl 

^nyn 


IIC/jvll < |v|i, VvGV/i. 


(5.60) 
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In view of Theorem 14.61 we can obtain the super convergence of the re¬ 
covered displacement gradients by following the same routine as in 

the proof of Theorem 4.2 in [53]. 

Theorem 5.1. Let u G and u/i G V/i 6e the displacement 

solutions of the problems (12.2p and (13.141) . respeetively. Under Diagonal 
condition (MCI) and Neighboring condition (MC2), the gradient reeov- 
ery is supereonvergent in the sense that 

||Vu - GhUhW < /i^+"(||u ||3 -h ||cj1|i) -h h'^Wcrh- (5.61) 


5.2 Recovery of stresses 

From the super convergence of the recovered displacement gradients 
in Theorem 15.11 we can easily derive the following super convergence of 
the recovered stresses (G/iU/j -|- for the stress tensor 

cr = Ce(u): 

||cr - GlohW < ||C|| (/i^+"(||u||3 -(- ||(t||i) -h hfWah) ■ (5.62) 

this estimate is not uniform with respect to 


However, due to the factor 
the Lame constant A. 

In what follows we shall construct a uniform recovered-type stress ap¬ 
proximation by following the idea of |46] . 

Denoting 

M/i := {u G L^(D) : D = o Fr-G span{l, for all FT G 7/i}, 

we introduce a recovered-type operator 

Rh : L^{D) ^ Mh 

as follows. For any i) G L^(D), we first define Rhi) at all vertices of 7h, then 
obtain Rhfj G on the whole domain by interpolation using the nodal 
shape functions of the piecewise isoparametric bilinear interpolation. 

For any interior vertex Zi{xi,yi) of 7h, we assume its patch uji, defined in 
(I5.56|> . consists of N elements, Ki,K 2 , • • • , Kj^, with N > 3. To define Rhtp 
at Zi we introduce the space W := span{l,x,y} and let (jff satisfy 


J{^f) = m^nJ{w), J{w) 


(5.63) 


Assume = ai + a 2 X -|- a^y and denote 


A,: = 


1 , 


'Ki 


Ik, 


Ik, 


2/ , 
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b : = 




then, from (|5.63l) . the constant vector a = ( 01 , 02 , 0 : 3 )^ is determined by 


A^Aa = A^h. (5.64) 

Thus it follows 

4>t = i'i-,x,y){A'^A) ^ A'^h. (5.65) 

We hence define 

Rh^iZo) = clA{Zo). (5.66) 

We next define Rh'4^ at any vertex Zf, E 50. Let Z^ be shared by m 
{m > 1) patches, e.g. a;i,a; 2 ,-- - ,iOm, which are corresponding m interior 
vertices Zi, Z 2 , • • • , Z^, then we can define 


- Hi 

Rhi^{Z,):=-J24>t{Zb), (5.67) 

i=\ 


where cl)f is given by (15.6511 . 

As a result, for any given stress finite element function r = (rn, T 22 , ti 2 )^ E 
S/i, we define the stress recovery R^t E M| with 

RhT := {RhTii, RhT22, RhTu)'^■ (5.68) 

Remark 5.2. We can show that A in (I5.64p is invertible for sufficiently 
small h. Since N > 3, it suffices to show rank{A) = 3. In fact, in view of 
(j3.3ll and (|3.1ip it holds 


A, = 


K 


Jl) *^0 


\K, 


+ A 


'1 


+4^4^), 


"0 


AT, 


+ 


'1 




II 

bi'‘lK,\)+0(hj+p, 

which implies 

( \Ki\ 

4^\Ki\ 

bS'lKil \ 

A = 

\K2\ 

4"\K2\ 



\ \Kn\ 

a4\KN\ 



+ 0(/l2+“), 


(5.69) 


where \Kj\ = 0{h‘^) is the area of Kj C wq. Recalling that Zq is an interior 
vertex of Th and ) is the center of the element Kj (1 < j < N,N > 

3 ), we easily have the fact that there exist at least three center points which 
are not lying on a same line. Thus, it holds rank(A) = 3 for sufficiently 
small h. 
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By the definition of Rh, we can derive Lemmas 15. 3115.41 
Lemma 5.3. The operator Rh '■ L‘^{Q) — Mh is bounded in norm with 

\\Rh^P\\<\m, VV^gL2(L!). (5.70) 

In addition, under Diagonal condition (MCI) it holds 

u - Run < e H^{n). (5.71) 

Proof. We first prove (|5.70n . Let V be the set of all vertices of Th- For 
^lJ € L^(n) and Zi{xi, pi) G V, let (ff € VF be the solution of the minimization 
problem (I5.63p . From (|5.66l) - (|5.67p we have 

WRhf^f « 

Zi&V 

Recalling that ||^||oo ^ ^ h~^ (cf. Remark l5.2p and b = 

(Iki /i ^2 ■ ■ ■ ’ fKpf (|5-b5p we easily obtain 

4>t{Zif = \{l,Xi,yi) {A^A) ^ A'^hf 

N 

i=i 

which, together with ^5.721) . leads to the desired conclusion. 

By noticing that the operator Rh preserves linear polynomials on each 
patch Wj, namely Rhi) = ip for ip G W, the desired estimate (15.7111 follows 
from the Bramble-Hilbert lemma and Diagonal condition (MCI). □ 

Lemma 5.4. For a G Ti,let G B/j be defined as in (14.2211 . Then it holds 

Rhcr = Rhcr^. (5.73) 

Proof. In light of (j5.66ll - (j5.67p . it suffices to show cp^^'- = cp^u. By the relation 
(j4.24ll it holds fj^ cru = for j = 1,2,-•• ,N. Then the conclusion 

follows from the minimization problem (15.6311 . □ 

Theorem 5.5. Let (a, u) G fl ^ i'^hi w) € 

B/iXV/i be the solutions of the problems (12.211 and p3.14p . respectively. Then, 
under Diagonal condition (MCI) and Neighboring condition (MC2), 
the following superconvergent result holds: 

Ik - RhC^hW < /i^’^"lklli + ^^Iklb- (5.74) 

Proof. By Lemma 15.41 it holds 

a - RhCfh = {<y - Rh<y) + Rh{(y^ - crh)- 

Then the desired superconvergence (I5.74p follows from Lemma (5 .3 1 and The¬ 
orem 14.61 □ 
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5.3 A Posteriori Error Estimates 

Denote e'^ := ||Vu — Vu/i||, e®" := ||cj — ah\\- Recall that Gh^h and 
Rh<^h are the recovered displacement gradients and the recovered stresses, 
respectively. In what follows we shall use the a posteriori estimators 

= \\GhUh-Vuh\\, = \\Rh(^h - (^h\\ 

to estimate the errors 

Theorem 5.6. Assume that Th satisfy Diagonal condition (MCI) and 
Neighboring condition (MC2). Let {a, a) € 

and (cr/i,u/i) G x Yh be the solutions of the problems (12.2h and (13.141) . 
respeetively. Then it holds 

- ||Vu - GhUhW <e^<g'^ + ||Vu - GhUh\\, (5.75) 

- ||(T - RhCThW < + ||cj - RhcrnW- (5.76) 

Moreover, if the solution ((T/i,u/i) is such that ||Vu — Vu/j|| > h and ||(t — 

'^h\\ ^ h, then the reeovery type a posterior error estimators are 

asymptotically exact in the sense 

rf/e'^ = 1 + g^/e^ = 1 + 0(/i““^“’i>). (5.77) 

Proof. The inequalities (|5.75p - (j5.76D follow from triangular inequality di¬ 
rectly, and the estimates (j5.77D follow from (|5.75p - (l5.76|l . Theorem 15.II and 
Theorem 15.51 

□ 


6 Numerical Experiments 

In this section we compute two test problems, Examples Ih.lllR^ to ver¬ 
ify our results of super convergence and a posterior error estimation for the 
PS hybrid stress finite element method. The examples are both plane strain 
problems with pure displacement boundary conditions, where the Lame pa¬ 
rameters g, A are given by 

E ^ Ev 

2(l + z/)’ “ (l + j/)(l-2z/)’ 

with 0 < < 0.5 the Poisson ratio and E the Young’s modulus. We set 

E = 1500. In all the computation we use 4x4 Gaussian quadrature. Notice 
that 2x2 Gaussian quadrature is accurate for computing the stiffness matrix 
of the PS hybrid stress FEM. All the fine meshes are obtained by bisection 
scheme. We compute the following relative errors for the displacement and 
stress approximation: 

iU ._ - Vwll 

|u|i |u|i 


9'^ := 


|W - u |i 
|u|i 


__ \Uh - U|i 

^ lull ' 
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Fh - 0- 


e := 


\\cr-ah\\ 7]^ \\Rh(^h-crh\ 

-, V ■= — “- 


||cr|| ||cr|| ||cr|| ||CJ| 

Example 6.1. The domain Q. = [0,1] x [0,1], the body force 

f=F'7r^( sin(vr?/) \ 

y — sin( 7 rx) cos( 7 ry) J ’ 

and the exact solution (u, a) is given by 

(1 + u) cos(7rx) sin(7ry) — 2(1 — v‘^)xy 


u = 


— (1 + v) sin(7rx) cos(7ry) + (1 — n‘^)x‘^ + v{\ + v){y‘^ — 1) 
—vr sin( 7 ra:) sin( 7 ry) — 2y 0 


a = E 


0 vr sin(7rx) sin(7r7/) 

The initial mesh is shown in Figure [67R and numerical results are listed in 
Table\^ 

Example 6.2. The domain Q. = [0,10] x [—1,1], the body force 

f = F-rr'^ ( sin(7r?/) \ 

y — sin(7rx) cos(7ry) J ’ 

and the exact solution is given by 

(1 + u) cos(7rx) sin(7r7/) — 2(1 — v‘^)xy 


u = 


— (1 + v) sin(7rx) cos(7ry) + (1 — i''^)x‘^ + v{\ + n){y‘^ — 1) 
—vr sin(7ra:) sm(7ry) — 2y 0 


a = E 


0 vr sin(7rx) sin(7r7/) 

The initial mesh is shown in Figure [UTSl and numerical results are listed in 
Table\^ 

We note that the refinement by bisection means that Diagonal con¬ 
dition (MCI) is satisfied with a = 1. From Tables [T][2] we can draw the 
following conclusions. 


• 0“ and 6 ^ are of second order convergence, uniformly with respect to 
A. These are conformable to the uniform superconvergence results in 
Theorem 14.61 

• e“ and fj^, as well as e'^ and f]^, are of first order convergence, uni¬ 
formly with respect to A. In particular, fj'^ and fj^ are asymptotically 
exact, which means the a posteriori estimators and are asymp¬ 
totically exact. All these are conformable to the a posterior estimates 
in Theorem 15.61 
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Figure 6.2: 2x2 irregular mesh for Example 16.11 
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Figure 6.3: 5x1 irregular mesh for Example 16.21 

References 

[1] G. Acosta and R. Duran. Error estimate for Qi isoparametric elements 
satisfying a weak angle condition. SIAM J. Numer. Anal., 38: 1073- 
1088, 2001. 

[2] D.N. Arnold, D. Boffi and R. Falk. Approximation of quadrilateral finite 
elements. Math. Comp. 71: 909-922, 2002. 

[3] I. Babuska, J. Oden and J. Lee. Mixed-hybird finite element approxima¬ 
tions of second-oder elliptic boundary-value problems. Comput. Meth¬ 
ods Appl. Mech. Engrg., 11(2): 175-206, 1977. 

[4] I. Babuska and T. Strouboulis. The Finite Element Method and its 
reliability. Oxford University Press, London, 2001. 

[5] 1. Babuska and M. Suri. On locking and robustness in the finite element 
method. SIAM J. Numer. Anal., 29(5): 1261-1293, 1992. 


21 







Table 1: The results of PS element on irregular meshes: Example 16.11 


V 

Error 

8x8 

16 X 16 

32 X 32 

64 X 64 

128 X 128 

Order 

0.3 


0.0051 

0.0013 

0.0003 

0.0001 

0.0000 

1.98 



0.1114 

0.0556 

0.0278 

0.0139 

0.0069 

1.00 


fjU 

0.1216 

0.0573 

0.0280 

0.0139 

0.0069 

1.03 



0.0138 

0.0034 

0.0009 

0.0002 

0.0001 

2.00 



0.0953 

0.0475 

0.0237 

0.0119 

0.0059 

1.00 



0.1059 

0.0491 

0.0240 

0.0119 

0.0059 

1.04 

0.49 

6'^ 

0.0054 

0.0014 

0.0004 

0.0001 

0.0000 

1.98 



0.1143 

0.0569 

0.0284 

0.0142 

0.0071 

1.00 


fjU 

0.1240 

0.0586 

0.0287 

0.0142 

0.0071 

1.03 



0.0153 

0.0038 

0.0010 

0.0002 

0.0001 

2.00 



0.1182 

0.0593 

0.0297 

0.0148 

0.0074 

1.00 



0.1294 

0.0609 

0.0299 

0.0149 

0.0074 

1.03 

0.4999 

6^ 

0.0057 

0.0014 

0.0004 

0.0001 

0.0000 

1.99 



0.1144 

0.0570 

0.0285 

0.0142 

0.0071 

1.00 


fjU 

0.1241 

0.0587 

0.0287 

0.0143 

0.0071 

1.03 



0.0155 

0.0039 

0.0010 

0.0002 

0.0001 

2.00 



0.1203 

0.0604 

0.0302 

0.0151 

0.0076 

1.00 


rf 

0.1315 

0.0620 

0.0304 

0.0151 

0.0076 

1.03 


[6] R.E. Bank and J. Xu. Asymptotically exact a posteriori error estima¬ 
tors, part I: grids with superconvergence. SIAM J. Numer. Anal., 41(6); 
2294-2312, 2003. 

[7] R. E. Bank and J. Xu. Asymptotically exact a posteriori error estima¬ 
tors, Part II: General unstructured grids. SIAM J. Numer. Anal., 41(6): 
2313-2332, 2003. 

[8] J. Brandts and M. Kfizek. Gradient super convergence on uniform sim- 
plicial partitions of polytopes. IMA Journal of Numerical Analysis, 23: 
489-505, 2003. 

[9] F. Brezzi and M. Fortin. Mixed and Hybird Finite Element Methods. 
Springer-Verlag, New York, 1991. 

[10] G.M. Chen. Structure Theory of Superconvergence of Finite Elements. 
Hunan Science Press (in Chinese), 2001. 

[11] L. Chen. Superconvergence of tetrahedral linear finite elements. Int. J. 
Numer. Anal. Model, 3: 273-282, 2006. 


22 
















5x1 


10 X 2 


Figure 6.4: Irregular meshes 


Table 2: The results of PS element on irregular meshes: Example 16.21 
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